Analysis of Measles data in Germany
From Ecological inference for infectious disease data, with application to vaccination strategies
21 June, 2019
1 Introduction and notes
- This file contains figures and analyses presented in ``Ecological inference for infectious disease data, with application to vaccination strategies“, along with additional analyses that are not discussed in detail in the primary manuscript.
- There are three MCMC’s run through Stan in this analysis, each may take approximately 30 minutes to compile.
2 Descriptive tables and figures
2.1 Measles cases and estimated vaccination coverage
Number of measles cases and estimated vaccination coverage for the 16 German states from 2005-2007. MMR1 (MMR2) is the estimated vaccination coverage for at least 1 (2) MMR vaccine, estimated from the school entry examinations. Note this partially reproduces Table 1 from Herzog et al 2011.
State Population Total.Cases MMR1 MMR2
1 Baden-Wuerttemberg (BW) 10,738,753 162 90.0% 75.6%
2 Bavaria (BY) 12,492,658 606 88.7% 73.2%
3 Berlin (BE) 3,404,037 104 90.0% 80.2%
4 Brandenburg (BB) 2,547,772 18 93.9% 86.9%
5 Bremen (HB) 663,979 4 88.4% 71.9%
6 Hamburg (HH) 1,754,182 29 90.0% 80.5%
7 Hesse (HE) 6,075,359 336 91.2% 78.1%
8 Mecklenburg-Western Pomerania (MV) 1,693,754 4 93.6% 88.0%
9 Lower Saxony (NI) 7,982,685 144 91.2% 78.0%
10 North Rhine-Westphalia (NW) 18,028,745 2,036 89.7% 76.9%
11 Rhineland-Palatinate (RP) 4,052,860 85 90.8% 77.3%
12 Saarland (SL) 1,043,167 0 91.0% 81.8%
13 Saxony (SN) 4,249,774 18 94.3% 82.4%
14 Saxony-Anhalt (ST) 2,441,787 12 94.1% 86.5%
15 Schleswig-Holstein (SH) 2,834,254 89 89.9% 79.3%
16 Thuringia (TH) 2,311,140 8 94.8% 85.9%2.2 Summary figures
2.3 Maps of the estimated proportion of individuals with at least one and two MMR vaccines.
2.4 Observed bi-weekly incidence of reported cases:
3 Fit hierarchical model in Stan
3.1 Get prior
priorch <- function(x,q1,q2,p1,p2){
(p1-pbeta(q1,x[1],x[2]))^2 + (p2-pbeta(q2,x[1],x[2]))^2 }
# p1 and p2 define the upper and lower levels of the interval
p1 <- 0.05; p2 <- 0.95
# q1 and q2 define the bounds of the interval
q1 <- 0.6; q2 <- 0.95
## together, we interpret these values to say:
## we want 90% (p2-p1) of the mass to be between q1 and q2
## solve and get the parameters for a beta distribution 90% of mass between 0.6 and 0.95
opt <- optim(par=c(1,1),fn=priorch,q1=q1,q2=q2,p1=p1,p2=p2)
opt$par[1] 10.004052 2.4490063.2 Fit model and assess converged
measlesdat=list(I=16, # number of areas
Ts=78, # number of time points
Nobs=16*78, # total number of observations
sin1=sin3, cos1=cos3, # inputs for seasonal terms
y=t(Y), # observed cases dim I x Ts
logNi=log(pop/sum(pop)), # log pop fraction
X=X1, # vaccination coverage
N=sum(pop), # total pop
a=c(10, 2.5) # parameters specifying prior on phi
)
pmt=proc.time()
fit <- stan(
model_code = mod_ENre,
chains=4, iter=4e3,
control=list(adapt_delta=0.99, max_treedepth = 15),
seed=47,
data = measlesdat)In file included from C:/Users/lfisher/Documents/R/win-library/3.5/BH/include/boost/config.hpp:39:0,
from C:/Users/lfisher/Documents/R/win-library/3.5/BH/include/boost/math/tools/config.hpp:13,
from C:/Users/lfisher/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core/var.hpp:7,
from C:/Users/lfisher/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
from C:/Users/lfisher/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core.hpp:12,
from C:/Users/lfisher/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/mat.hpp:4,
from C:/Users/lfisher/Documents/R/win-library/3.5/StanHeaders/include/stan/math.hpp:4,
from C:/Users/lfisher/Documents/R/win-library/3.5/StanHeaders/include/src/stan/model/model_header.hpp:4,
from file2cec7ae77eb9.cpp:8:
C:/Users/lfisher/Documents/R/win-library/3.5/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
# define BOOST_NO_CXX11_RVALUE_REFERENCES
^
<command-line>:0:0: note: this is the location of the previous definition(proc.time()-pmt)/60 user system elapsed
0.04483333 0.03866667 17.28033333 print(fit, pars=c("alpha_ar", "phi", 'alpha_en', 'gamma', 'delta', 'sigma_ar', 'sigma_en', "lp__"), probs=c(.025,.5,.975)) # Inference for Stan model: c4c1ff19f358c6c9ee8fb5649733aef4.
4 chains, each with iter=4000; warmup=2000; thin=1;
post-warmup draws per chain=2000, total post-warmup draws=8000.
mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
alpha_ar 0.85 0.01 0.50 -0.26 0.91 1.66 1166 1
phi 0.89 0.00 0.09 0.66 0.92 0.99 2508 1
alpha_en 3.47 0.01 0.43 2.54 3.52 4.15 1956 1
gamma 0.71 0.00 0.08 0.55 0.71 0.86 8000 1
delta -0.20 0.00 0.08 -0.36 -0.20 -0.04 8000 1
sigma_ar 0.74 0.01 0.35 0.28 0.66 1.61 685 1
sigma_en 0.55 0.00 0.17 0.28 0.52 0.95 1963 1
lp__ 8516.14 0.17 5.65 8504.09 8516.47 8526.34 1153 1
Samples were drawn using NUTS(diag_e) at Mon Apr 29 16:18:41 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).print(fit, pars=c("alpha_ar", 're_ar', 'sigma_ar', "lp__"), probs=c(.025,.5,.975)) # Inference for Stan model: c4c1ff19f358c6c9ee8fb5649733aef4.
4 chains, each with iter=4000; warmup=2000; thin=1;
post-warmup draws per chain=2000, total post-warmup draws=8000.
mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
alpha_ar 0.85 0.01 0.50 -0.26 0.91 1.66 1166 1
re_ar[1] 0.43 0.01 0.31 -0.07 0.38 1.18 795 1
re_ar[2] 0.55 0.01 0.31 0.08 0.50 1.29 679 1
re_ar[3] 0.30 0.01 0.32 -0.22 0.26 1.07 755 1
re_ar[4] 0.30 0.01 0.43 -0.49 0.28 1.21 1838 1
re_ar[5] -0.62 0.01 0.70 -2.28 -0.51 0.46 3029 1
re_ar[6] -0.80 0.01 0.55 -2.14 -0.71 0.03 3129 1
re_ar[7] 0.64 0.01 0.30 0.19 0.60 1.38 709 1
re_ar[8] -0.26 0.01 0.65 -1.78 -0.20 0.86 5237 1
re_ar[9] 0.20 0.01 0.32 -0.32 0.16 0.95 1017 1
re_ar[10] 0.76 0.01 0.30 0.32 0.71 1.49 671 1
re_ar[11] 0.20 0.01 0.34 -0.37 0.16 0.96 1094 1
re_ar[12] 0.02 0.01 0.87 -1.69 0.01 1.80 6090 1
re_ar[13] -0.44 0.01 0.54 -1.74 -0.37 0.48 4084 1
re_ar[14] -0.93 0.01 0.70 -2.60 -0.81 0.10 2232 1
re_ar[15] 0.37 0.01 0.32 -0.15 0.33 1.12 959 1
re_ar[16] -0.74 0.01 0.71 -2.46 -0.62 0.33 2522 1
sigma_ar 0.74 0.01 0.35 0.28 0.66 1.61 685 1
lp__ 8516.14 0.17 5.65 8504.09 8516.47 8526.34 1153 1
Samples were drawn using NUTS(diag_e) at Mon Apr 29 16:18:41 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).print(fit, pars=c("alpha_en", 're_en', 'sigma_en', "lp__"), probs=c(.025,.5,.975)) # Inference for Stan model: c4c1ff19f358c6c9ee8fb5649733aef4.
4 chains, each with iter=4000; warmup=2000; thin=1;
post-warmup draws per chain=2000, total post-warmup draws=8000.
mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
alpha_en 3.47 0.01 0.43 2.54 3.52 4.15 1956 1
re_en[1] -0.26 0.00 0.22 -0.70 -0.25 0.17 3013 1
re_en[2] 0.34 0.00 0.22 -0.08 0.33 0.80 2408 1
re_en[3] 0.58 0.00 0.25 0.11 0.57 1.09 2588 1
re_en[4] -0.32 0.00 0.33 -1.02 -0.31 0.30 8000 1
re_en[5] -0.10 0.00 0.36 -0.84 -0.09 0.60 8000 1
re_en[6] 0.61 0.01 0.27 0.09 0.60 1.16 2505 1
re_en[7] 0.00 0.00 0.25 -0.51 0.00 0.50 3438 1
re_en[8] -0.56 0.01 0.39 -1.40 -0.54 0.13 4267 1
re_en[9] 0.37 0.00 0.22 -0.06 0.36 0.81 2547 1
re_en[10] -0.01 0.00 0.21 -0.42 -0.02 0.41 2514 1
re_en[11] 0.45 0.00 0.25 -0.04 0.45 0.95 2989 1
re_en[12] -0.85 0.01 0.50 -2.01 -0.79 -0.04 3283 1
re_en[13] -0.41 0.00 0.29 -1.01 -0.40 0.13 3819 1
re_en[14] -0.04 0.00 0.29 -0.62 -0.03 0.52 4432 1
re_en[15] 0.50 0.01 0.26 0.00 0.49 1.02 2688 1
re_en[16] -0.24 0.00 0.32 -0.90 -0.23 0.35 8000 1
sigma_en 0.55 0.00 0.17 0.28 0.52 0.95 1963 1
lp__ 8516.14 0.17 5.65 8504.09 8516.47 8526.34 1153 1
Samples were drawn using NUTS(diag_e) at Mon Apr 29 16:18:41 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).traceplot(fit, pars=c("alpha_ar", "phi", 'alpha_en', 'gamma', 'delta','sigma_ar', 'sigma_en', "lp__")) #
pairs(fit, pars=c("alpha_ar", "phi", 'alpha_en', 'gamma', 'delta', 'sigma_ar', 'sigma_en', "lp__")) #
3.3 Parameter estimates
50% 2.5% 97.5%
alpha_ar 0.911 -0.256 1.656
phi 0.917 0.659 0.989
alpha_en 3.525 2.543 4.155
gamma 0.708 0.551 0.860
delta -0.200 -0.356 -0.043
sigma_ar 0.663 0.277 1.607
sigma_en 0.524 0.279 0.955
R0 2.488 0.774 5.2373.4 Examine model fits
3.5 Estimates of random effects
3.6 Examine posterior distributions
3.7 Plot Pearson residuals
4 Additional analyses
4.1 Comparison of MLE vs Bayesian analysis
4.1.1 MLE estimates
loglik_eco=function(param, y, xi, n){
v=expit(param[2])
l=0
nwks=dim(y)[2]
ni=dim(y)[1]
for(i in 1:ni){
for(t in 2:nwks){
lam=(exp(param[1])*y[i, t-1] + (n[i]/sum(n))*exp(param[3] + param[4]*sin(2*pi*t/26) + param[5]*cos(2*pi*t/26)) )
l=l-dpois(y[i, t], lambda=(1-v*xi[i])*lam, log=T)
}
}
return(l)
} parameter Est SE
1 alpha_ar 2.02928493 0.07239919
2 phi 0.98867494 0.68640739
3 alpha_en 4.06273036 0.10100270
4 gamma 0.69281062 0.08418688
5 delta -0.07449201 0.081201284.1.2 Estimates from fixed effects Stan model
measlesdat2=measlesdat
measlesdat2$a=c(1, 1) # put a flat prior on phi
pmt=proc.time()
fit_fe <- stan(
model_code = mod_fe,
chains=4, iter=4e3, # warmup = 1e3, # save_dso = T, # warmup=7e3, # refresh=F, # #
control=list(adapt_delta=0.99, max_treedepth = 15),
data = measlesdat2)In file included from C:/Users/lfisher/Documents/R/win-library/3.5/BH/include/boost/config.hpp:39:0,
from C:/Users/lfisher/Documents/R/win-library/3.5/BH/include/boost/math/tools/config.hpp:13,
from C:/Users/lfisher/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core/var.hpp:7,
from C:/Users/lfisher/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core/gevv_vvv_vari.hpp:5,
from C:/Users/lfisher/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/core.hpp:12,
from C:/Users/lfisher/Documents/R/win-library/3.5/StanHeaders/include/stan/math/rev/mat.hpp:4,
from C:/Users/lfisher/Documents/R/win-library/3.5/StanHeaders/include/stan/math.hpp:4,
from C:/Users/lfisher/Documents/R/win-library/3.5/StanHeaders/include/src/stan/model/model_header.hpp:4,
from file29ac349f2e0d.cpp:8:
C:/Users/lfisher/Documents/R/win-library/3.5/BH/include/boost/config/compiler/gcc.hpp:186:0: warning: "BOOST_NO_CXX11_RVALUE_REFERENCES" redefined
# define BOOST_NO_CXX11_RVALUE_REFERENCES
^
<command-line>:0:0: note: this is the location of the previous definition(proc.time()-pmt)/60 user system elapsed
0.06216667 0.04533333 3.48383333 print(fit_fe, pars=c("alpha_ar", "phi", 'alpha_en', 'gamma', 'delta'), probs=c(.025,.5,.975), digits_summary = 3) # Inference for Stan model: 8a16ff8514700ac19fc4f3422c38a24a.
4 chains, each with iter=4000; warmup=2000; thin=1;
post-warmup draws per chain=2000, total post-warmup draws=8000.
mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
alpha_ar 2.066 0.002 0.057 1.918 2.080 2.141 1380 1.003
phi 0.993 0.000 0.007 0.975 0.996 1.000 1307 1.003
alpha_en 4.120 0.002 0.087 3.922 4.126 4.273 2068 1.002
gamma 0.698 0.001 0.084 0.534 0.699 0.863 5049 1.001
delta -0.074 0.001 0.080 -0.232 -0.074 0.082 4526 1.000
Samples were drawn using NUTS(diag_e) at Thu May 02 14:29:44 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).4.2 Sensitivity analysis - measles analysis with non-informative priors
measlesdat=list(I=16, # number of areas
Ts=78, # number of time points
Nobs=16*78, # total number of observations
sin1=sin3, cos1=cos3, # inputs for seasonal terms
y=t(Y), # observed cases dim I x Ts
logNi=log(pop/sum(pop)), # log pop fraction
X=X1, # vaccination coverage
N=sum(pop), # total pop
a=c(1, 1) # parameters specifying prior on phi
)
pmt=proc.time()
fit_flat <- stan(
model_code = mod_ENre,
chains=4, iter=4e3, # warmup = 1e3, # save_dso = T, # warmup=7e3, # refresh=F, # #
control=list(adapt_delta=0.99, max_treedepth = 15),
data = measlesdat)
(proc.time()-pmt)/60 user system elapsed
0.01650000 0.02766667 21.40100000 print(fit_flat, pars=c("alpha_ar", "phi", 'alpha_en', 'gamma', 'delta', 'sigma_ar', 'sigma_en', "lp__"), probs=c(.025,.5,.975)) # Inference for Stan model: c4c1ff19f358c6c9ee8fb5649733aef4.
4 chains, each with iter=4000; warmup=2000; thin=1;
post-warmup draws per chain=2000, total post-warmup draws=8000.
mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
alpha_ar 0.97 0.03 0.77 -0.91 1.24 1.85 687 1
phi 0.87 0.01 0.22 0.15 0.97 1.00 892 1
alpha_en 3.58 0.02 0.70 1.83 3.85 4.32 871 1
gamma 0.71 0.00 0.08 0.55 0.71 0.87 7386 1
delta -0.20 0.00 0.08 -0.36 -0.20 -0.03 6306 1
sigma_ar 0.69 0.01 0.33 0.24 0.62 1.53 772 1
sigma_en 0.51 0.00 0.18 0.24 0.49 0.95 1567 1
lp__ 8521.75 0.21 6.10 8509.33 8522.09 8532.99 841 1
Samples were drawn using NUTS(diag_e) at Mon Apr 29 16:44:34 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).print(fit_flat, pars=c("alpha_ar", 're_ar', 'sigma_ar', "lp__"), probs=c(.025,.5,.975)) # Inference for Stan model: c4c1ff19f358c6c9ee8fb5649733aef4.
4 chains, each with iter=4000; warmup=2000; thin=1;
post-warmup draws per chain=2000, total post-warmup draws=8000.
mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
alpha_ar 0.97 0.03 0.77 -0.91 1.24 1.85 687 1
re_ar[1] 0.38 0.01 0.31 -0.11 0.34 1.11 792 1
re_ar[2] 0.50 0.01 0.31 0.03 0.45 1.22 707 1
re_ar[3] 0.27 0.01 0.32 -0.24 0.23 1.01 849 1
re_ar[4] 0.30 0.01 0.42 -0.51 0.28 1.15 1992 1
re_ar[5] -0.58 0.01 0.66 -2.25 -0.47 0.45 3012 1
re_ar[6] -0.75 0.01 0.51 -1.91 -0.68 0.03 3393 1
re_ar[7] 0.61 0.01 0.29 0.17 0.57 1.28 793 1
re_ar[8] -0.23 0.01 0.61 -1.62 -0.16 0.80 4131 1
re_ar[9] 0.18 0.01 0.31 -0.33 0.14 0.88 959 1
re_ar[10] 0.72 0.01 0.29 0.28 0.67 1.41 728 1
re_ar[11] 0.17 0.01 0.32 -0.38 0.14 0.90 1108 1
re_ar[12] -0.01 0.01 0.75 -1.56 -0.01 1.50 6256 1
re_ar[13] -0.39 0.01 0.52 -1.59 -0.32 0.47 4406 1
re_ar[14] -0.84 0.02 0.68 -2.47 -0.73 0.15 1849 1
re_ar[15] 0.33 0.01 0.32 -0.19 0.29 1.07 894 1
re_ar[16] -0.66 0.01 0.66 -2.24 -0.55 0.32 1964 1
sigma_ar 0.69 0.01 0.33 0.24 0.62 1.53 772 1
lp__ 8521.75 0.21 6.10 8509.33 8522.09 8532.99 841 1
Samples were drawn using NUTS(diag_e) at Mon Apr 29 16:44:34 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).print(fit_flat, pars=c("alpha_en", 're_en', 'sigma_en', "lp__"), probs=c(.025,.5,.975)) # Inference for Stan model: c4c1ff19f358c6c9ee8fb5649733aef4.
4 chains, each with iter=4000; warmup=2000; thin=1;
post-warmup draws per chain=2000, total post-warmup draws=8000.
mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
alpha_en 3.58 0.02 0.70 1.83 3.85 4.32 871 1
re_en[1] -0.27 0.00 0.23 -0.72 -0.27 0.19 2620 1
re_en[2] 0.30 0.01 0.23 -0.13 0.28 0.78 1654 1
re_en[3] 0.54 0.01 0.26 0.07 0.53 1.05 2166 1
re_en[4] -0.27 0.01 0.33 -0.97 -0.25 0.33 3655 1
re_en[5] -0.11 0.00 0.36 -0.86 -0.10 0.58 5456 1
re_en[6] 0.56 0.01 0.28 0.02 0.55 1.12 2166 1
re_en[7] 0.00 0.00 0.25 -0.49 0.00 0.48 3239 1
re_en[8] -0.51 0.01 0.39 -1.39 -0.47 0.15 3284 1
re_en[9] 0.35 0.00 0.22 -0.07 0.35 0.80 2379 1
re_en[10] -0.03 0.00 0.21 -0.42 -0.04 0.39 2050 1
re_en[11] 0.42 0.01 0.26 -0.06 0.41 0.94 2378 1
re_en[12] -0.79 0.01 0.50 -1.97 -0.72 -0.01 3016 1
re_en[13] -0.36 0.01 0.30 -0.98 -0.34 0.18 3299 1
re_en[14] -0.01 0.01 0.29 -0.60 0.00 0.53 3066 1
re_en[15] 0.46 0.01 0.26 -0.03 0.45 1.01 2159 1
re_en[16] -0.20 0.01 0.32 -0.88 -0.18 0.40 3619 1
sigma_en 0.51 0.00 0.18 0.24 0.49 0.95 1567 1
lp__ 8521.75 0.21 6.10 8509.33 8522.09 8532.99 841 1
Samples were drawn using NUTS(diag_e) at Mon Apr 29 16:44:34 2019.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).traceplot(fit_flat, pars=c("alpha_ar", "phi", 'alpha_en', 'gamma', 'delta','sigma_ar', 'sigma_en', "lp__")) #
pairs(fit_flat, pars=c("alpha_ar", "phi", 'alpha_en', 'gamma', 'delta', 'sigma_ar', 'sigma_en', "lp__")) #
4.2.1 Examine model fits
4.2.2 Estimates of random effects
4.2.3 Examine posterior distributions
5 Stan Models
Below, we print the code to define the Stan models fit in this document
mod_ENre='
data{
int<lower=1> I; // number of areas
int<lower=1> Ts; // number of times
int<lower=1> Nobs; // number of total observations
int y[I, Ts] ; // IxTs matrix of observations
real logNi[I]; // log population proportions
vector[I] X; // vector of Xi vacc cover
vector[Ts] sin1; // vector of sin
vector[Ts] cos1; // vector of cos
vector[2] a; //vector of parameters for the beta prior on phi
}
parameters{
real alpha_ar; // AR component - intercept
real alpha_en; // EN component - intercept
real gamma; // sin component
real delta; // cosine component
real<lower=0, upper=1> phi; // vaccine effect
vector[I] re_en; // EN random effects
vector[I] re_ar; // AR random effects
real<lower=0> sigma_ar;
real<lower=0> sigma_en;
}
model{
for(i in 1:I){
re_ar[i]~normal(0, sigma_ar);
re_en[i]~normal(0, sigma_en);
for(t in 2:Ts){
real muit;
muit = exp(alpha_ar+re_ar[i])*y[i, t-1] + exp(logNi[i] + alpha_en + re_en[i] + gamma*sin1[t] + delta*cos1[t]);
y[i, t] ~ poisson((1-phi*X[i])*muit);
}
}
//priors
alpha_ar ~ normal(0, 5);
alpha_en ~ normal(0, 5);
gamma ~ normal(0, 10);
delta ~ normal(0, 10);
phi ~ beta(a[1], a[2]);
}
generated quantities{
vector[Nobs-I] log_lik;
vector[I] ar;
vector[I] en;
vector[I] r;
for(i in 1:I){
r[i] = exp(alpha_ar + re_ar[i])*(1-phi*X[i]);
ar[i] = exp(alpha_ar + re_ar[i]);
en[i] = exp(alpha_en + re_en[i]);
for(t in 2:Ts){
real muit;
muit=exp(alpha_ar + re_ar[i])*y[i, t-1] + exp(logNi[i] + alpha_en + re_en[i] + gamma*sin1[t] + delta*cos1[t]);
log_lik[(i-1)*(Ts-1)+t-1 ] = poisson_lpmf(y[i, t] | (1-phi*X[i])*muit);
}
}
}
'
### fixed effects model - to the MLE
mod_fe='
data{
int<lower=1> I; // number of areas
int<lower=1> Ts; // number of times
int<lower=1> Nobs; // number of total observations
int y[I, Ts] ; // IxTs matrix of observations
real logNi[I]; // log population sizes
vector[I] X; // vector of Xi vacc cover
vector[Ts] sin1; // vector of sin
vector[Ts] cos1; // vector of cos
vector[2] a; //vector of parameters for the beta prior on phi
}
parameters{
real alpha_ar; // AR component - intercept
real alpha_en; // EN component - intercept
real gamma; // sin component
real delta; // cosine component
real<lower=0, upper=1> phi; // vaccine effect
}
model{
for(i in 1:I){
for(t in 2:Ts){
real muit;
muit = exp(alpha_ar)*y[i, t-1] + exp(logNi[i] + alpha_en + gamma*sin1[t] + delta*cos1[t]);
y[i, t] ~ poisson((1-phi*X[i])*muit);
}
}
//priors
alpha_ar ~ normal(0, 5);
alpha_en ~ normal(0, 5);
gamma ~ normal(0, 10);
delta ~ normal(0, 10);
phi ~ beta(a[1], a[2]);
}
generated quantities{
vector[Nobs-I] log_lik;
vector[I] ri;
for(i in 1:I){
ri[i] = exp(alpha_ar)*(1-phi*X[i]);
for(t in 2:Ts){
real muit;
muit=exp(alpha_ar)*y[i, t-1] + exp(logNi[i] + alpha_en + gamma*sin1[t] + delta*cos1[t]);
log_lik[(i-1)*(Ts-1)+t-1 ] = poisson_lpmf(y[i, t] | (1-phi*X[i])*muit);
}
}
}
'